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Abstract 

This paper concerns the numerical study for the generalized Rosenau-Kawahara-RLW equation obtained by coupling 
the generalized Rosenau-RLW equation and the generalized Rosenau-Kawahara equation. We first derive the energy 
conservation law of the equation, and then develop a three-level linearly implicit difference scheme for solving the 
equation. We prove that the proposed scheme is energy-conserved, unconditionally stable and second-order accurate 
both in time and space variables. Finally, numerical experiments are carried out to confirm the energy conservation, 
the convergence rates of the scheme and effectiveness for long-time simulation. 
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1. Introduction 

The nonlinear wave is one of the most important scientific research areas. During the past several decades, many 
scientists developed different mathematical models to explain the wave behavior, such as the KdV equation, the RLW 
equation, the Rosenau equation, and many others. In the following, we give a short review of these important wave 
models. 

The well known KdV equation 

M/ + Mx.tA + = 0, (1) 

was first introduced by Boussinesq [1] in 1877 and rediscovered by Diederik Korteweg and Gustav de Vries [2] in 
1895. Since then, there are a lot of studies on this equation and its variational form. Here we just mention some of 
the recent work. Kudryashov [3] reviewed the travelling wave solutions for the KdV and the KdV-Burgers equations 
proposed by [4], Biswas [5] studied the solitary wave solution for KdV equation with power-law nonlinearity and time- 
dependent coefficients, while Wang et al. [6] investigated the solitons, shock waves for the potential KdV equation. In 
addition to the theoretical studies, readers can refer to [7, 8] for the numerical simulations of the KdV equation and 
the generalized KdV equation. 

The regularized long-wave (RLW) equation (also known as Benjamin-Bona-Mahony equation) 

Ut + Ux + UUx — Uxxt — 0, (2) 

was first proposed as a model for small-amplitude long wave of water in a channel by Peregrine [9, 10]. The regular¬ 
ized long-wave (RLW) equation and its different variational forms were well studied both theoretically and numeri¬ 
cally in the literature. Readers can refer to [11-17] for numerical studies and [18-20] for theoretical studies. 

When study the compact discrete systems, the well-known KdV equation can not describe the wave-wave and 
wave-wall interactions. To overcome the shortcoming of KdV equation, Rosenau proposed the following so-called 
Rosenau equation [21, 22]: 

^xxxxt ^x UUjl — 0. (3) 
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The existence and uniqueness of the solution for the Rosenau equation were theoretically proved by [23]. Besides the 
theoretical analysis, numerical studies of equation (3) also exist in the literature, see [24-28] and references therein. 

For further consideration of nonlinear waves, the term -u „t is included in the Rosenau equation. The resulting 
equation is usually called the Rosenau-RLW equation [29-31]; 


Ut - + Mxx.vxf + Mx + MMx = 0. (4) 

The above equation was further extended into the generalized Rosenau-RLW equation [32-34]: 

~ ^^xxt ^^xxxxt + Ux + M Ux — 0, (5) 

where m > 1 is a positive integer. The initial boundary value problems for the Rosenau-RLW equation and the 
generalized Rosenau-RLW equation have been well studied numerically in the past years [29-34]. For example. 
Pan et al. [29] developed a three-level linearly implicit conservative scheme for the usual Rosenau-RLW equation (4), 
where the method achieves second-order accuracy both in time and space variables. Zuo et al. [32] proposed a second- 
order Crank-Nicolson finite difference scheme for the generalized Rosenau-RLW equation (5), however, the method is 
a nonlinear scheme, iterations are needed at each time step. Later on. Pan et al. [33] used a three-level linearly implicit 
conservative scheme for the same generalized Rosenau-RLW equation (5), where the method is a linear scheme and 
is also second-order convergent both in time and space variables. 

On the other hand, to consider another behavior of nonlinear waves, the viscous term Mxxx needs to be included in 
the Rosenau equation (3). The resulting equation is usually called the Rosenau-KdV equation [35-39]: 

Ut + Uxxxxt + Uxxx + Ux + UUx = 0. (6) 

For numerical investigations, Hu et al. [39] proposed a second-order linear conservative finite difference method for 
the Rosenau-KdV equation. However, numerical methods for the initial-boundary value problem of the Rosenau-KdV 
equation have not been studied widely. 

By coupling the above Rosenau-RLW equation and Rosenau-KdV equation, one can obtain the following Rosenau- 
KdV-RLW equation [40-44], 

Ut - Uxxt + Uxxxxt + Uxxx + Ux + UUx = 0. (7) 

For numerical investigation, Wongsaijai et al. [40] proposed a three level implicit conservative finite difference method 
for the above Rosenau-KdV-RLW equation. For theoretical studies, the solitary wave, shock waves, conservation laws 
as well as the asymptotic behavior for the Rosenau-KdV-RLW equation with power law nonlinearity were studied 
by [41^4], where the power law nonlinearity means the the last term in the left-hand side of equation (7) is replaced 
by a general nonlinear term (m^)x and p is any positive integer. 

In addition, the following Kawahara equation 

Ut + Ux -t" UUx "t" Uxxx ~ Uxxxxx — O5 (8) 

arose in the theory of shallow water waves with surface tension [45]. Equation (8) is called the modified Kawahara 
equation if the third nonlinear term in the left-hand side is replaced by u^Ux- There is a wide range of literature on the 
numerical investigations and theoretical studies for the usual Kawahara equation and the modified Kawahara equation. 
For theoretical aspects, some periodic and solitary wave solutions for both the Kawahara equation and the modified 
Kawahara equation are provided in [46^8]. In addition to the theoretical studies, readers can refer to [49-51] for the 
numerical studies of the Kawahara equation and the modified Kawahara equation. 

As one more step consideration of the nonlinear wave, Zuo [38] obtained the Rosenau-Kawahara equation by 
adding another viscous term —Uxxxxx to the Rosenau-KdV equation (6), and studied the solitary solution and periodic 
solution of the Rosenau-Kawahara equation. The Rosenau-Kawahara equation is given as follows [38]: 

Ut + Ux + UUx + Uxxx + Uxxxxt “ Uxxxxx = 0. (9) 

For numerical study, Hu et al. [52] proposed a two level nonlinear Crank-Nicolson scheme and another three-level 
implicit linear conservative finite difference scheme for the Rosenau-Kawahara equation, both methods are proved to 
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be second-order convergent. And Biswas [53] investigated the solitary solution and the two invariance of following 
equation, 

Uj -r CLllx + Ux + CUxxx ^^XXXXt ^^XXXXX — Os (10) 

where a, b, c, v, a, A are real constants, m is a positive integer, which indicates the power law nonlinearity. Equation 

( 1 0) is referred as the generalized Rosenau-Kawahara equation [52] . As far as we know, besides the analytical solitary 
solution obtained by [53], there is no numerical study for this generalized Rosenau-Kawahara equation. 

Besides all the literature mentioned in the above, there are also a lot of related references which used both an¬ 
alytical and numerical methods to study different nonlinear models. For example, Triki et al. [54] discussed the 
solitary solution for the Gear-Grimshaw model, Bhrawy et al. [55] provided solitons, cnoidal waves and snoisal 
waves for the Whitham-Broer-Kaup system, while Shokri et al. [58] numerically solved the two-dimensional complex 
Ginzburg-Landau equation by using the meshless radial basis functions method, and Dehghan et al. [59] applied a 
semi-analytical method for solving the Rosenau-Hyman equation. Readers can refer to [54—57] for more discussions 
on analytically finding solitons and other wave solutions for different nonlinear models, and refer to [58-60] for more 
discussions on numerical studies. 

In this paper, we will study the following generalized Rosenau-Kawahara-RLW equation: 

U, + aUx + bu"’Ux + CUxxx - aUxxt + AUxxxxt - Vllxxxxx -0, Xl< X< Xr, 0 <t <T. (11) 

which is obtained by coupling the generalized Rosenau-RLW equation (5) and the generalized Rosenau-Kawahara 
equation (10). Analytical solitary solutions of the above equation (11) were recently obtained by [61] through the 
sech ansatz method. 

In this work, we shall consider the numerical simulation of above equation (11) with initial condition 

u(0,x) - Uq(x), Xl < X < Xr, (12) 

and boundary conditions 

u{xi, f) - u{Xr, t) - 0, Ux(xi, t) - Ux{Xr, t) - 0, Uxx(Xl, t) - UxxiXr, t) - 0, 0 < t < T, (13) 

where m is a positive integer, a, b, c, v are all real constants, while a, A are set to be positive constants, xi is a large 
negative number and Xr is a large positive number. In this paper, we only discuss the solitary solution of equation 

(11) . By solitary wave assumptions, the solitary solution and its derivatives have the following asymptotic values: 
M —> 0 as X —» + 00 , and ^ 0 as x ^ +oo, for « > 1. Thus, the boundary conditions (13) are meaningful for the 
solitary solution of equation (11). And we assume that the wave peak is initially located at x = 0, and x;, Xr, which are 
large numbers, are used to assure that the solitary wave peak always locates inside the domain [x/, Xr] during the time 
interval [0, T], Similar set-ups are used in [39, 40, 52]. 

Note that equation (11) reduces to the generalized Rosenau-RLW equation when c - Q,v - 0. Specially, it reduces 
to the usual Rosenau-RLW equation (4) when a = m - a = A - l,b - 2,c = Q,v - 0. And equation (11) reduces 
to the generalized Rosenau-Kawahara equation when a = 0. Specially, it reduces to the usual Rosenau-Kawahara 
equation when a -m-c-A-v - \,b - 2,a - 0. 

For numerically solving differential equations, the total accuracy of a particular method is affected not only the 
order of accuracy of the method, but also other factors. The conservative property of the method is another factor 
that has the same or possibly even more impact on results. For example, one successful and active research is to con¬ 
struct structure-preserving schemes (or called symplectic schemes) for the ODE systems (see [62] and the references 
therein). Better solutions can be expected from numerical schemes which have effective conservative properties rather 
than the ones which have nonconservative properties [63]. And Li et al. [64] even pointed out that in some areas, the 
ability to preserve some invariant properties of the original differential equation is a criterion to judge the success of 
a numerical simulation. 

For the initial boundary value problem (ll)-(l3), we will show that it satisfies a fundamental energy conservation 
law in the next section. In addition, the wave equation is highly nonlinear due to the third nonlinear term in the 
left-hand side of the equation (11). When considering the finite difference scheme for the equation (11), the usual 
Crank-Nicolson scheme will lead to a nonlinear scheme with heavy computation while other standard linearized 
discretizations for the nonlinear term, e.g., one step Newton’s method or a second-order extrapolation method, will 
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loss the the energy conservative property. An ideal scheme should have relative less computational cost, can preserve 
energy, be unconditionally stable and maintain second-order accuracy. 

In this paper, a three-level linearly implicit finite difference method for the initial value problem (11)-(13) will be 
presented. The fundamental energy conservation is preserved by the presented numerical scheme. The existence and 
uniqueness of the numerical solution are also proved. Moreover, numerical analysis shows that the method is second- 
order convergent both in time and space variables, and the method is unconditionally stable. Numerical results confirm 
well with the theoretical results. 

The rest of the paper is organized as follows: Section 2 shows the energy conservation of the initial boundary value 
problem (11)-(13). Section 3 gives the detail description of the three-level linearly implicit finite difference method, 
the proof for the discrete conservative property, the existence and uniqueness as well as the convergence and stability 
of the numerical solution. Numerical results are shown in section 4. Conclusions are provided in the final section. 


2. Conservative property 

Equations (11)-(13) satisfy the following energy conservative property. 

Theorem 1. Suppose uq e Cq[x;, x^], then the solution of (ll)-(l3) satisfies the energy conservation law: 


E{t) 


-f"- 


(x, t) + airfix, t) + Au\fix, t)dx 




u^(x, 0) -I- au\{x, 0) -H Au\fix, Q)dx - E{Q), 


(14) 


for any t e [0, T], where Cg [x/, Xr] is the set of functions which are seventh order continuous differentiable in the 
interval [x;, Xr] and have compact supports inside (x/, xfi. 

Proof. Multiplying (11) by 2m and integrating over the interval [x/, x^], one get 


d 2 T' r 

— I udx + a I 2uu^dx + b ( 

dt J^, J^, J^, 

; I 2uu_xxtdx + A I 

Xj Xi 


2uuxdx+ b 2u'”^Uu)xdx + c 2uuxxxdx 


a 2uuxxtdx + A 2uUxxxxtdx - V '^UUxxxxxdx - 0. 


f 

*J X[ 

r 


(15) 


Using the integration by parts, one can easily obtain 

r' r' 1/2 2 ^ 

I UUxdx — I udu — — lu (Xr, t) — U (xi, f)) = 0, 

kJxi *J Xi ^ 


f u”'^^{u)xdx- f 

Jx, Jx, 


u”'^‘du = 


1 


m + 2 


\xi 


-U 


= 0 , 


/" 

X] 

r 

Xj 


UUxxxdx — MM, 


Ullxxtdx - MM, 


iv, r' r' r 

I - I Uxxdu - — I UxxUxdx - - I 

Jxj J Xi *J Xj 

Id 2, id r 

2 dt Jx, 2 dt Jx, 


1 


UxdUx — -M?!'' = 0, 

JC 2 Ax, 


u dx. 


^ ^^Exxxfdx — ^ udUxxxl — ^^^xxxt\^^ ^ ^xtlxxxtdx — UxUxxt\„^^ "i" ^ ^xxd^ — 2 ^ ^^xxd^-> 

r xj - 1 -'^ - r xt -- r ^ _ 1 2 1 -. _. 

I ^^ExxxxUX — UUxxxx\,„^ I ^xxxxtlxttX — UxUxxx\„^^ "t" I ^xxxtlxxU:!: — x^^xxlx, — 

A Xi d X] A Xi ^ 

where the boundary conditions (13) are used. 

Thus, only the first, fifth and sixth term in the left hand side of (15) are nonzero, all other terms are zero. This 
yields. 


dtJx, " 


^ -I- aul + Aul,.dx — 0. 


(16) 
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Therefore, 


E{t) - I i?‘{x, f) + au\{x, t) + Au\j^x, f)dx — I i?{x, 0) + au^ix, 0) + Au^^(x, 0)dx - E(0), 

*Jxi *Jxi 


(17) 


for any t 6 [0, T], This completes the proof. 

3. Numerical method 

3.1. Numerical scheme 

In this section, we give a complete description of our numerical method for the initial value problem (11)-(13). We 
first describe our solution domain and its grid. The computational domain is defined as Q = {(x, f)|x; < x < x^, 0 < 
t < T}, which is covered by a uniform grid Q/, = {(x,-, f„)|x, = x/ + ih, f„ - nr, i = 0, • • • , M, n = 0, • • • , A^), with 
spacing h - t We denote t/" as the numerical approximation of m ( x ,', t„) and Z° - {U - (t/,)|t/_i - Uq - 
U\ - Um-i - Um = Um+i =0,1 = -1,0,1,-'' ,M- \,M,M+ 1). For convenience, the difference operators, inner 
product and norms are defined as follows; 


t7" = 


+ U’}- 


f/«+l _ 


iui), = 


2t 

M-l 




W U"-U" , 

h ^ ^ h 


{Ulh = F") = /z 2 t/fFf, II U" ||2= {U\ V), II U" iu= , max _ \U"\. 


1=1 




where t/" = (t/"j, • • • , U'h^i) is the numerical solution at time t„ - m. 

The essential of our scheme is that the third term in the left-hand side of (11) is rewritten and discretized as 

bu"'u, = -\(iru, + ^ -±-[(uT(0f)sc + ((uTum, 

m + 2 m + 2 


(18) 


and which is a second-order approximation around (x, - xi + ih, = nr). And all other terms of left-hand side of (11) 
are discretizated by using the standard second-order central difference method around (x, - xi + ih, t„ - m). 

The detailed numerical scheme is as follows: 


and 


{UD, + a(U"h + - :,[(UT(Uih + ((UTUih] + c(Ufh.n - a{Ul),„ 

m + 2 

+ T(t/f),,a, - v(t7")xxm = 0, 2<i<M-2, \<n<N-\. 
t/f = uoixi), 0<i<M, 

Ui {Ui), = (t/^)i = 0, {Ui),-, = (t/^),^ = 0, 0 < ; < A. 


(19) 

( 20 ) 
( 21 ) 


Obviously, the above conditions (21) will give U\ - = 0 on two internal points and = 0 on two 

fictitious points, for any 0 < j < N. Thus, t/l 6 Z®, for any 0 < j < N. In addition, we can see that (19) is a three-level 
linearly implicit scheme and the coefficient matrix of the linear system (19) is banded. Therefore, the resulting linear 
system of equations can be solved efficiently using a direct linear solver, such as the LU decomposition method. 

Since the scheme is a three-level method, to start the computation, we need to give the method for computation of 
t/'. The t/* is computed through the following Crank-Nicolson scheme: 


T 

/■rrl 


+ a 


'U] + t/f' 

b 

'Uj + t/f' 

m 

iuUu^\ 

'Uj + U°' 

m 

'Uj + t/f' 

\ 

2 J 

m + 2 

X 

2 J 


[ 2 J, 

2 J 


i 2 J 

h. 


U\ - t/f' 


U} + t/° 


+ A 


U - U\ 




ft/.' H- t/°) 


= 0 , 


( 22 ) 
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where it is a nonlinear scheme and is second-order accurate both in time and space variables. 

We point out that this discretization (18) for the third term in the left-hand side of (11) is specially designed so 
that the numerical scheme (19)-(21) can guarantee the energy conservation (see the proof of theorem 2 in the next 
subsection). Readers can refer to references [29, 40, 52] for the similar kind of treatments for the nonlinear term uujc in 
their equations to achieve conservative schemes. In addition, our scheme (19)-(22) is only nonlinear at first time step 
when computing U\ and all other steps are linear. Thus, the computational cost is relatively cheap. In the following 
several sections of the paper, we will show that our scheme is second-order accurate both in time and space variables, 
can preserve the energy identity (14), and is unconditionally stable. 

The following lemmas are well known results, which are essential for existence, uniqueness, convergence, and 
stability of our numerical solution. Lemma 1 and 2 can be verified through direct computation, lemma 3 is the dis¬ 
crete form of the Sobolev’s embedding theorem, which can be found in lemma 1, page 110 of [65]. In the rest part 
of the paper, unless otherwise indicated, C is the notation referring to a general positive constant, which may have 
different values in different contexts. 


Lemma 1. For any two mesh functions U,V € Z®, one have 

(Ui, V) = -{U, Vst), (U„ V) = -{U, V,), {U,-„ V) = -{U„ V,), 


Furthermore, 

iu, U,,,,) =11 U,, ||2 . 

Lemma 2. For any mesh function U 6 Zj^’, one have 

(t/i, U) = 0, (t/,a, U) = 0, (t/„m, u) = 0. 

Lemma 3. (Discrete Sobolev’s inequality (Lemma 1, page 110 of [65]) For any mesh function U 6 Zj^’, one have 

II U iu< c II U, II . 


3.2. Discrete conservation 

Theorem 2. Suppose mq g C'^ix^Xr], then the solution of finite difference scheme (19)-(21) satisfies || U” ||oo< C, 
II U" ||oo< C, for any 0 < n < N, and also satisfies the following discrete energy conserx’ation: 


A II t/""' r + II t/" IP II t/f ‘ IP + II U’] IP , II t/"L 11 ^ + II t/L 

E = -r-H cr---h A --- 


rO ||2 


II t/' IP + II IP II t/i IP + II IP , II UE II- + II t/L „ ^ 

"- 2 - “- 2 - ^ - 2 -^ ^ ’ 

for any 0 < n < N — 1, where E” is the discrete energy at time t — (n + i)T. 


(23) 


Proof. Taking the conditions - t/g = uj = = ^m+i = 0 (0 ^ j ^ N) into account, and after 

computing the inner product of equation (19) with U", i. e., vve have 


rrf7+l ||2 _ II Tjn-\ ||2 l 

.-H a(ui, u") + ^ ((t/")'"(t7")jE -H ((U'TU'fs, 0”) 


4t ' ^ ' m + 2 

+ c(t7L,, t7") - a(U:,„ t7") + A(U:,,,„ U") - U") = 0. 


(24) 

(25) 


By using lemma 2, we get 


(U"„ U’f = 0, (t/L., t7") = 0, U") = 0. 


(26) 
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Moreover, 


and 


([(t/"r (t7"). + [{U'T t7")J, t7") 

, M-\ 

= T Z [iU'iriUli - UU) + - {Ul^'Ul^ U’l 

^ (=1 

, M-l , M-1 

= ^ I] liu’iru’/^.u’i + (ui,rui,u"] --J] \(uij"ufuu + (u-ruiui^ 


= 0 . 


u:p II" - II II' 

4t 


(U" - U") - 

xxxxt’ ^ / 


, 11^-liar'll' 


4t 


( 27 ) 


(28) 


where boundary conditions (20) and (21) are used. 
Thus, 


II [/«+i ||2 _ II t/«-' ||2 +a II t/r' II' -a II t/r' II' +A II t/”7‘ -A || ||2= 0, 

for any 1 < n < - 1. This is equivalent to 

II ||2 + II t/« ||2 II t/f 1 ||2 + II U: ||2 II t/«+' ||2 + II t/«- ||2 

-;;-+ o: -r-H A ---= 

2 2 2 

II U" ||2 + II t/«-' ||2 II t/« ||2 + II t/”-l ||2 II f/«. ||2 + II [/-I ||2 

-2--2--2-’ 

for any 1 < n < - 1. This further yields 

E"^E^, for any 1 <11 < N - I, 

which is actually the energy conservation law (23). 

u'+ifi 

Multiplying (22) both sides by ' ^ ' and using the similar techniques as above, one can obtain 
II t/“ ||2 +a II U° II' +A II U% ||2=|| t/‘ ||2 +a || Ul ||' +A || t/], ||2 . 

Thus, (31) can be rewritten as 


(29) 


(30) 


(31) 


(32) 


E" =11 ||2 +a II t/“ II' +A II ||2 . (33) 

Since mq 6 Cq[x/, and the initial condition (20) are used in the numerical method, the right-side of (33) is 
bounded. By assumptions, a, A are positive constants, therefore, 

II U" ||< C, II 11^ C, for any 0 < n < M (34) 

By using lemma 3, we have || f/" ||oo< C. 

In addition, through direct computation, one can verify that 

WU'iJHlU’i.W, forany0<n<M (35) 


Thus, 


II t/l- ||< C, for any 0<n<N. (36) 

Again by using lemma 3, we have || U" ||oo< C. This completes the proof. 
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3.3. Existence and uniqueness 

Theorem 3. The finite difference scheme (19)-(21) has a unique solution. 


Proof. To prove the theorem, we proceed by the mathematical induction. Suppose 1/', ■ ■ ■ , t/"(l < n < N - 1) 
are solved uniquely, we now consider the equation (19) for Assume that are two solutions of (19) 

and let jj jg g^gy jq verify that satisfies the following equation; 




-nf-l 




2{m + 2) 




+ = o. 


(37) 


Taking the inner product of (37) with W"^', we have 




-nf-l ||2 


+ ^ II 

2t 


-nf-l ||2 


II ii'= 0, 


(38) 


where 


W"^^) = -(ITf', ITf'), =11 Ip, {W’lfl,, = 0, 

(ITf',lT"+‘) =0, ([(t/T(H^''^‘).+((t/TH^''^‘y =0, (39) 

are used. The first five identities of (39) are directly from lemma 1 and 2, and the last one can be obtained as follows: 


([(t/")”' + {(u'T 

. M-\ 

= 2 Z - ^i-x) +- iu’urwti] wr^ 


/=i 

M-l 


M-\ 


= 5 Z ^ Xi Yuirwt^wf^ + iu’urwt^wr^] 


1=1 


= 0 . 


(40) 


From (38) and the definition of the || ■ ||-norm, one can see that (38) has only a trivial solution. Thus, (19) 
determines t/"^' uniquely. This completes the proof. 


3.4. Convergence and stability 

Let u(x,t) be the solution of problem (11)-(13), t/" be the solution of the numerical schemes (19)-(21), and 
u" = u{xi, tn), e'l = u" - U”, then the truncation error of the scheme (19)-(21) can be obtained as follows; 

b 


dl = {ef), + aie’D, + ^ [{uyfuf), + («)'" u'}), - (U^f (U-), - ((t/f)'" U’})^ 


m + 2 

^ ^xxx ~ 4(^i )xxxxt - y (C; 


(41) 


where e" = ' , 2 < i < M - 2 and 1 < « < A - 1. 

Since all terms in (19) are the second-order approximations of the corresponding terms in left-hand side of (11) 
around (x, - xi + ih, t„ = nr), by Taylor expansion, it can be easily obtained that r" = 0(t^ + h^) if /i, t —> 0 and 
u(x, t) 6 where the set of functions which are seventh order continuous differentiable in space and third 

order continuous differentiable in time. This following lemma is a well known result. 
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Lemma 4. (Discrete Gronwall’s inequality). Suppose that w(k) and p(k) are nonnegative functions while p(k) is a 
non-decreasing function. If 

k-\ 

w(k) < p(k) + Cr ^ w(l), V k, 

1=0 


then 


w(k) < V k. 

Theorem 4. Suppose uq 6 Cq[xi, Xr], and u{x, t) 6 C^'^, then the numerical solution t/" of the finite difference scheme 
(19)-(21) converges to the solution of the problem (11)-(13) in the sense of \\ ■ ||oo, and the convergence rate is 
0 (t^ + h^X i.e., 

II u" - U" |U< C(t^ + h\ for any 2<n<N. 


Proof. Taking the inner product of (41) with 2e", we have 

II ||2 - II e"-' ||2 +a(|| ef* Ip - || ll') + ^(11 11^ - II iP) 

= 2t [(r", 2e”) - a (e"„ 2e") - c ((e\,,, 2e") + v , 2e") - (Q, 2e") - (R, 2e")], 


where 


Q 


[{u'T (jn, - iuT (t/"),], R = [{(u'T - {(UT t/")J. 


(42) 


(43) 


(44) 


m + 2 

By using lemma 2, we obtain 

(4,2e''}^0, ((e'%,,,2el^0, ((e'Vm, 2e") = 0. (45) 

From the notations introduced at the beginning of section 3, for any 0 < n < N, we have the following inequality 


«z i z z K^r * m 

j=i j=i j=i 


1=1 1=1 
Thus, for any 0 < n < N, we have 


h =11 e' 


n ||2 


(46) 


In addition, we have 

KG, 2e«)| 


2bh 


m + 2 
2bh 




1=1 

M-l 


m + b 
2bh 


^ [(«p'”(spi - (up'"(t7;),. + (up^Kt/pi - (t/p'-Kt/p^] e] 


y=i 

M-l 


Z <";)”<':)<^+44 Z lb:)" - wpl wp**: 


M-l 


m + 2 
2bh 


1=1 

M-l 




1=1 

M-l tm-l 


m + 2 


1=1 


1=1 ik=0 




<ChY,[\{e4k\ + \e"^]r;\ 

1=1 

< C(|| ef' 11^ + II ep‘ ir + II e"' ‘ IF + II e" IF + II e" . 

< C(|| ef' ||2 + II ep' ||2 + II ||2 + II e’' |p + || e"-' ||p 


( 47 ) 
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where theorem 2, Cauchy-Schwarz inequality and inequality (46) are used. 
Similarly, we have 




M-l 


m + 2 

2bh 


M-l 


2bh 


m + 2 

2bh 

m + 2 


j=i 

M-l 

z 

>1 

M-l 



'm-l 

-7 


^ k=0 } Si 



7=1 




k=0 


ie’’)s 




7=1 

< C(|| ef' Ip + II 4 -' ||2 + II ||2 + II ||2 + II ||2) 


„«-l ,|2 ^ II g«+l ||2 II ||2 ^ II g«-l ||2x 


< C(|| r + II 4 

where theorem 2, Cauchy-Schwarz inequality and inequality (46) are used again. 

Furthermore, we have 

(r«, 2^") < 2 II r" |||| e" ||<|| r" ||2 + || g" ||2<|| r" Ip +1 (|| ||2 + || e'-' ||2), 

where Cauchy-Schwarz inequality are used. 

Substituting (47)-(49) into (43), we get 


e«+l ||2_ II g«-l ||2 


+a(|| ef* ||2 - II e'r' Ih + ^(11 II' - II ||') 

< Ct(|| II' + II e" II' + II e"-' II' + II ef ‘ ||' + || ||') + 2t || r" ||' . 

Since a, A are positive constants, it is easy to check that 

II II'- II e”-' II' +a(|| ef' II' - II t-r' II') + 4(11 II" - II II") 


< Ct(|| ef* II' + II e" 


^n+l ||2 


e"ll' + lk"-' ||') + 2T||r" II' 


< C'rdl g 


.n+l Ii2 


+2 II e" II' + II 


-1 ii2 


+cr||er' II' +2a||e^||^) + 2T||r" II 


n || 2 \ 


where C' - max(-, C) and C is the positive constant in the above inequality. 

Replacing C' in the above inequality by the general positive constant notation C, we have 


^,.rl ||2_ II ^«-l ||2 


+a(|| ef' II' - II er' ||') + 4(|| ||' - || e";' ||') 

< Ct(|| II' +2 II e" II' + II e"-' ||' +a || ef' ||' +2a || ||') + 2 t || r" ||' 

< Ct(|| II' +2 II e" ||' + || e"-' ||' +a || ef* ||' +2a || ||' +a || ||' 

+ A II II' +24 II e”, ||' +4 || ||') + 2t || r" ||' . 


Let 


(48) 


(49) 


(50) 


( 51 ) 


D" =11 e" II' +a II e:: II' +4 II eL ||2 + || g-i ||2 +£, || 


+4 II ||2 
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then (51) can be rewritten as follows: 


(£)"+' - D") < Ct(D'’+' + D") + 2t II r" H^, 


which is equivalent to 

(1 - Ct)(D'’+' - D") < ICtD" + 2t II r" . (52) 

If r, which is sufficiently small, satisfies t < ^ (C is the positive constant in the inequality (52)), then 1 - Cr > 0 
and (52) gives 

£)"+i _ D" < -2£—tD" + II r" ||2 

1 - Ct 1 - Ct 

< 3CtD" + 3t II r" 

< C"(tD" + t II r" ih, 

where C" - max(3C, 3) and we have used < 3 since t < ^. 

Replacing C" in the above inequality by the general positive constant notation C, we have 


£)"+i _ £)« < CtD" + Ct II r" . 

(53) 

Summing (53) from 1 to « - 1, we get 


n-1 n—1 

D" < + CtYj Ik'll^ 

/=1 l=l 

(54) 

where 


«-i 

T y II r' ||k nr max || r' ||k T ■ C(t + h^f. 

^ KKn-l 

l=l 

(55) 

Since -Q and the Crank-Nicolson scheme (22) is used to compute we have 

simple analysis for the scheme (22). Therefore 

= 0{t^ + h^) followed by a 

n-l 

D" < 0{t^ + l?f + CtY^dK 

/-I 

(56) 

Using lemma 4, we obtain 


D" < Ofr^ + h^f. 

(57) 

Thus, 


Ik" ||< Oir^' + h^ 

(58) 

By using lemma 3, we have 


II k IU< 0{t^ + h\ 

(59) 

i.e., 


II m" - t/" iu< C(t2 + h^). 

(60) 


This completes the proof. 


Theorem 5. Suppose uq € Cq[xi, Xr], then the solution U" of the finite difference scheme (19)-(21 ) is unconditionally 
stable with the || ■ ||oo norm. 

The proof of this theorem is similar as the above theorem. 
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4. Numerical results 


In order to preform numerical accuracy test of the method proposed in this paper, one needs to know the exact 
solitary solutions beforehand. Exact solitary solutions of the equation (11) were previously obtained by [61] using the 
sech ansatz method. In this paper, we provide a different derivation through the sine-cosine method, which is listed in 
the appendix A. 

Example 1. We present the numerical results for the case m - 2, a - \,b - \,c — 2,a - l,A- l,v = 1. From 
[61] and appendix A, one can obtain the exact solitary solution as follows; 


3 V3^-5 VTO . 2 


u{x, 0 = T I - 

V5V37- 


-sech 


29 


VV37-5 

\ 

( 33-5V37) 

4 

V 

X r 

l 5V37-29 j 


and the initial condition is set as 


^ 3 V3^ - 5 VTO ^2 

£/(x, 0) = -— — seen 


a/5V37- 


29 


VV37-- 


( 61 ) 


(62) 


We first carry out the numerical convergence studies. For the spatial convergence, we set r = 0.005 as the fixed 
time step and use 4 different spatial meshes: h - 0.8,0.4,0.2,0.1, where t is sufficient small such that the temporal 
error is negligible comparing to the spatial error (Here the time step is t = 0.005 while the smallest spatial size is 
h - 0.1, thus !i^ » T~, therefore, the dominant errors are the spatial errors). The final time T is set to be 10, and 
xi - -40, Xr — 200. Table 1 gives the errors between numerical solutions and exact solutions. We can see that the 
error decreases when the spatial mesh is refined and the convergence rate is two. Thus, the method is second-order 
convergent in space, which is consistent with theoretical results in the above section. For the temporal convergence, 
we set h - 0.005 as the fixed spatial mesh and use 4 different temporal meshes; r = 0.8,0.4,0.2,0.1, where h is 
sufficient small such that the spatial error is negligible comparing to the temporal error (Here the smallest time step is 
T = 0.1 while the spatial size is h - 0.005, thus « t^, therefore, the dominant errors are the temporal errors). The 
final time T is set to be 10, and x; = -40, x^ = 200. Table 2 gives the errors between numerical solutions and exact 
solutions. Again, we can see that the error decreases when the temporal mesh is refined, and the convergence rate is 
also two. Thus, the method is second-order convergent in time, which is again consistent with theoretical results in 
the above section. 

In order to show that the numerical scheme has the energy conservative property (23), we carry out another 
computation, where T - 100, x/ = -40, x^- = 240,/z = 0.1, t = 0.1 are used. Table 3 gives the quantities of E" at 
several time stages. We can see that E" is conserved exactly (up to 8 decimals) during the time evolution of the solitary 
wave. 

Example 2. We present the numerical results for the case m - A, a - \,b - \,c — 2,a - 1,A- l,v = 1. From 
[61] and appendix A, one can obtain the exact solitary solution as follows; 


40( VT27- 10)2 ] 

1 

4 

sech 

-y/vn?- 10 

\ 

( 118-10VT27'1 

l3(ioVm- i09)J 

3 

\ 

r 10Vl27- 109'j 

/ 


and the initial condition is set as 


m(x, 0) 


' 40( Vl27- 10)2 Y 
,3(10 Vm- 109), 


sech 



(63) 


(64) 


Again, we carry out the spatial and temporal convergence. Table 4 and 5 give the errors between numerical 
solutions and exact solutions for spatial and temporal convergence, respectively. Once again, we can see that the 
method is second-order convergent both in time and space variables. 

Additionally, table 6 provides the quantities of £" with T - 100, x; = -40, x^ = 240, h - 0.1,t - O.l. Once again, 
we can see that E" is conserved and the method can be well used to study the solitary wave at long time. 
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5. Conclusions 


In this paper, a three-level linearly implicit finite difference method for the initial value problem of the generalized 
Rosenau-Kawahara-RLW equation is developed. The fundamental energy conservative property is preserved by the 
current numerical scheme. The existence and uniqueness of the numerical solution are also proved. The method is 
shown to be second-order convergent both in time and space variables, and it is unconditionally stable. Moreover, 
exact solitary solutions are derived through sine-cosine method which is used for the numerical tests. Numerical 
results confirm well with the theoretical results. 
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Table 1: Spatial mesh refinement analysis with t = 0.005, T = 10 for example 1. 


h 

II HI 

rate 

Ik Hoc 

rate 

0.8 

2.66e-l 

- 

1.032e-l 

- 

0.4 

6.650-2 

2.002 

2.570-2 

2.006 

0.2 

1.666-2 

1.997 

6.460-3 

1.992 

0.1 

4.209-3 

1.985 

1.631-3 

1.986 
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Table 2: Temporal mesh refinement analysis with h = 0.005, T = 10 for example 1. 


T 

II HI 

rate 

II dU 

rate 

0.8 

1.523 

- 

6.060e-l 

- 

0.4 

3.812e-l 

1.999 

1.540e-l 

1.976 

0.2 

9.711e-2 

1.973 

3.910e-2 

1.977 

0.1 

2.442e-2 

1.992 

9.820e-3 

1.993 
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Table 3: Invariant of for example 1. 


t 

E" 

0.05 

25.451405792697514 

19.95 

25.451405792693116 

39.95 

25.451405792447929 

59.95 

25.451405792214793 

79.95 

25.451405791920855 

99.95 

25.451405792207414 
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Table 4: Spatial mesh refinement analysis with t = 0.005, T = 10 for example 2. 


h 

II HI 

rate 

II dU 

rate 

0.8 

1.543e-l 

- 

5.839e-2 

- 

0.4 

3.790-2 

2.026 

1.446e-2 

2.013 

0.2 

9.440-3 

2.005 

3.599e-3 

2.007 

0.1 

2.366-3 

1.996 

9.01 le-4 

1.999 
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Table 5: Temporal mesh refinement analysis with h = 0.005, T = 10 for example 2. 


T 

II HI 

rate 

II dU 

rate 

0.8 

4.582-1 

- 

2.029e-l 

- 

0.4 

8.633e-2 

2.408 

3.843e-2 

2.401 

0.2 

2.124e-2 

2.023 

9.447e-3 

2.024 

0.1 

5.263e-3 

2.012 

2.330e-3 

2.019 
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Table 6: Invariant of for example 2. 


t 

E" 

0.05 

13.565665615099391 

19.95 

13.565665614771643 

39.95 

13.565665614965912 

59.95 

13.565665614937172 

79.95 

13.565665614960499 

99.95 

13.565665614998375 
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Appendix A. Exact solitary solutions 

The sine-cosine method uses the sine or cosine function as the wave function to seek the traveling wave solution of 
a time dependent partial differential equation, which has the advantage of reducing the nonlinear problem to a system 
of algebraic equations that can be easily solved by using a symbolic computation system such as Mathematica or 
Maple [3,38,40, 66 ], 

For equation (11), one can obtain the exact solitary solutions by using the sine-cosine method. Firstly, we seek the 
following travelling wave solution [3, 38, 40, 66 ]: 

u{x,t) = u(^, ^-x-vt, (A.l) 


where v' is referred as the wave velocity which is a constant to be determined later. 
Under the transformation of (A.l) and integrating once, (11) can be reduced into: 


(a - v)u + 


m + \ 


+ (va + c)u^^ - (Av + = 0. 


(A.2) 


Secondly, by using the sine-cosine method, the solutions of the above reduced ODE equation can be expressed in the 
following form [38, 40, 66 ]: 

' Acos''(B^), 

0 , 


fi(^) - 


otherwise, 


(A.3) 


or in the form 


m(^) = 


A sin''(B^), 

0 , 


otherwise. 


(A.4) 


where A, B, rj are parameters to be determined. 
Using (A.3), one have 


= AB-r]{ri - 1) cos''^2(Bf) - AB^rf cos^B^), 


(A.5) 


and 


= AB"^t](j] - l)(r] - 2 )(t] - 3) cos’’-'^(B^) - lAB'^rjiri - \){rf -Irj + l) cos’’^^(B^) + AB'^rf' cos''(B^). (A.6) 

Substituting (A.3), (A.5) and (A.6) into (A.2), one obtain 

-cos''''"+')(B^) - AB^rjir] - \)(ri - 2 )(ti - 3)(4v -H v) cos’’-'^(B^)+ 

m + I 

{AB~(va + c)ri(ri - 1) H- 2AB‘^ri{ri - IXif -2rj + 2)(Av + v)) cos’’^^{B^)+ 

[{a - v)A - AB\va + c)rf - AB^'rf'iAv + v)) cos''(Bf) = 0. (A.l) 

Balancing cos''*^'"^'^(B^) and cos''^"^(Bf), and setting each coefficients of cos^(B^) (y = 77 , 77 - 2, 77 - 4) to be zero, one 
can obtain a set of equations for 77 , v. A, B as follows: 

77(777 -I- 1) = 77 — 4, (A. 8 ) 

bA"'^^ 

-^ _ AB^rjir^ - \){tj - 2)(jj - 3)(4y + v) = 0, (A.9) 

777 -1-1 

AB^va + c)t](t] - 1 ) -h 2AB\(ri - 1)(t]^ -2t] + 2))(Av + v) ^ 0, (A. 10) 

(a — v)A - AB^(va + c)tj^ - AB‘^rf(Av H- v) = 0. (A. 11) 

From (A.8), one can find that 

77 = -i. (A.12) 

m 
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From (A. 10) and (A.l 1), one obtain that 




(Aa + v)(ri^ -2t] + 2) ± yj{Aa + v)^{Tf- - 2ti + 2)^ + (Ac - va)(aa + c)if-{ri - 2)^ 


(Ac - va)rf(ri - 2)^ 


Once B is obtained, one can get that 


V - - 


2vB‘^(ri- -2ri + 2) + c 
2ABHTf- -2ri + 2) + a 


A = 


m + 1 


-BV^-l)(/7-2)(77-3)(dy + y) 


The classification for the above solution is discussed as follows: 
Case 1: (Aa + v)^(ri^ — 2ri + 2)^ + (Ac — va)(aa + c)if'(ri — 2)^ > 0. 


( 1 ): 


(/la+v)(7/^-27/+2)+ (Aa+v'f {rf-~lrf+TA+{Ac~va){aa+c)if-(r]-2)^ 


{Ac-va)rfiri-iy- -^ obtain a periodic solution as follows: 

u(x, t) - u(^) - A cos'^(B(x — vt)). 


(A.13) 


(A. 14) 
(A. 15) 


(A. 16) 


where B is given by 


B = 


I (Aa + v)(t]^ — 2ri + 2) + ■\J(Aa + v)^(t]^ — 2ri + 2)^ + (Ac — va)(aa + c)rf-(ri — 2)^ 


(Ac - va)rf-(ri - 2)^ 


(A. 17) 


and rj, v,A are given by (A. 12), (A. 14) and (A. 15), respectively. 

\c-va){aa-¥c)j]-{ri-2y „ 

- < 0. V 

u(x, t) - u(^) - Asech''(Z?o( 2 i: - vt)). 


{Aa+v)(rr-2r]+2)+^J{Aa+v)-(TA-2rt+2f+(Ac~va)(aa+c)rf-(r]--2Y . 

(2): -— (Ac-va)rf{ri- 2 )'’- - obtain a solitary solution as follows: 


(A. 18) 


where Bq is given by 




I -(Aa + v)(rf- — 2ri + 2) — yj(Aa + v)~(rp- -2?] + 2)^ + (Ac - va)(aa + c)if-(ri - 2)^ 

(Ac - va)rp-(ri - 2)^ 


(A.19) 


and Tj, v,A are given by (A. 12), (A. 14) and (A. 15), respectively. 

(3) : - - — (Ac-va)ri'^{tj- 2 y- -^ Similar periodic solution as (A. 16) can be obtained. 

(4) : ^ ^ obtained. 

Case 2: (Aa + v)^(if — 2t] + 2)^ + (Ac — va)(aa + c)if'(rj - 2)^ < 0. 

In this case, the resulting A, B, v' are complex numbers with nonzero real parts. After substituting the expressions 
of A, B, V, T] into equation (A.3), u needs to take the real part or the imaginary part of the resulting equation. The detail 
is omitted here. 

Similarly, one can use the sine method to get another set of solutions. 
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